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1. Introduction. Algorithms for reordering sparse matrices play a vital role in our 
ability to perform many large-scale matrix. computations. Ordering algorithms such as 
minimum- degree and nested dissection have been developed for reducing fill in direct 
methods for solving sparse, symmetric positive definite systems of equations [7, 12, 27]. 
Various ordering algorithms for reducing the envelope (variable band or profile) of 
sparse matrices, such as the reverse Cuthill-McKee (RCM), Gibbs-Poole-Stockmeyer 
(GPS), and Gibbs-King (GK) algorithms, have also been designed [12, 15, 21]. Al- 
though envelope-reducing orderings were developed for use in envelope schemes for 
direct factorization, these orderings have been used in the past few years in several 
other applications. The RCM ordering has been found to be an effective preordering in 
computing incomplete factorization preconditioners for preconditioned conjugate gra- 
dients methods [6, 8]. Such orderings have also been used in parallel matrix- vector 
multiplication and tridiagonalization of sparse symmetric matrices. 

The wider applicability of envelope-reducing orderings justifies a fresh look at the 
algorithms currently available and the development of new algorithms. In this paper we 
present a new spectral algorithm for computing an envelope-reducing ordering of sparse, 
symmetric matrices. The ordering algorithm uses an eigenvector corresponding to the 
smallest positive eigenvalue of the discrete Laplacian matrix associated with the given 
symmetric matrix. (If the matrix is irreducible, or equivalently if its adjacency graph is 
connected, then this eigenvector corresponds to the second smallest eigenvalue. Hence 
we call this a second Laplacian eigenvector or Fiedler vector.) The ordering is computed 
by permuting the components of a second Laplacian eigenvector in nonincreasing (or 
nondecreasing) order. For large matrices, the eigenvector computation is performed by 
a ‘multilevel’ approach described in [3]. 

Earlier, we had used a second eigenvector of the Laplacian matrix for computing 
a spectral nested dissection ordering, and for partitioning computations on finite ele- 
ment meshes on a distributed-memory multiprocessor [29, 30, 31]. The eigenvector of 
the adjacency matrix corresponding to the largest eigenvalue has been used to find a 
pseudoperipheral node by Grimes et al. [17]. 

A companion paper [13] provides theoretical justification for the spectral envelope- 
reduction algorithm by considering a closely related problem called the 2-sum problem. 
(This problem is defined in the next section.) It is shown there that this problem 
can be formulated as a quadratic assignment problem involving the Laplacian matrix. 
Lower bounds for the 2-sum are obtained in terms of the smallest positive Laplacian 
eigenvalue. These bounds appear to be reasonably tight, and hence indicate how close 
the computed orderings are to the optimal orderings. Further, permuting the matrix 
in nonincreasing (or nondecreasing) order of the components of a second Laplacian 
eigenvector is shown to yield a feasible solution to the 2-sum problem that is closest to 
an infeasible solution for which the lower bound is attained. 

Fiedler [9, 11] studied the properties of the second Laplacian eigenvalue and a 
corresponding eigenvector and their relationship to the connectivity of a graph, and 
also observed [10] that the differences in the components of this eigenvector is an ap- 
proximate measure of the distance between the vertices. Juvan and Mohar [19] have 
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advocated the use of this eigenvector to compute bandwidth and p-sum reducing or- 
derings. Mohar and Poljak [25] have recently provided a comprehensive survey of the 
applications of Laplacian spectra to combinatorial problems. 

The spectral envelope-reduction algorithm has several features which set it apart 
from the earlier reordering algorithms such as the GPS, GK, or RCM algorithms [5, 12, 
15, 21]. These algorithms employ local-search in the adjacency graph of the matrix. All 
of them try to find a pseudo- diameter in the graph by generating a long level-structure 
by breadth-first-search beginning from a suitable vertex. These types of algorithms 
generally do not vectorize, and there is no obvious way to implement them in parallel. In 
contrast the new algorithm proposed here is based on the computation of an eigenvector 
of a special matrix, and hence involves standard floating point operations, such as 
matrix vector multiplications, dot products etc. The algorithms for these operations 
not only vectorize easily, but also can be implemented in parallel with little effort. 
(Parallel implementation of the basic spectral method, which uses the Lanczos algorithm 
to find eigenvectors, is straightforward. Parallel implementation of the ‘multilevel’ 
enhancements described in Section 3 is more difficult, but possible in principle.) The 
algorithm is also iterative in nature, in the same sense that SOR or the Lanczos methods 
are iterative. It allows a user to terminate the reordering process depending on a 
stopping criterion, thus permitting the user to make trade-offs in ordering time versus 
storage efficiency. 

Before we end this introduction, some comments axe in order about the applicabil- 
ity of the results to envelope factorization schemes. Frontal methods related to envelope 
or profile schemes are still the method of choice for solving large-scale systems of linear 
equations in many structural engineering applications, for example in the computa- 
tional structural mechanics testbed (CSM) at NASA Langley [20]. Implementations of 
these methods are also widely distributed in most of the finite element software pack- 
ages such as MSC/NASTRAN or ANSYS. Parallel algorithms for the actual numerical 
factorization of a matrix in envelope format have been investigated [28, 33]. 

Efficient implementations of sparse matrix algorithms [1, 2, 22, 32] on supercom- 
puters demonstrate that very high levels of performance are attainable with general 
sparse algorithms. Hence there are no good reasons to use envelope schemes for sparse 
matrix factorizations for the sake of performance alone. Furthermore, it has long been 
known that general sparse methods axe considerably more efficient with respect to stor- 
age [12]. Ashcraft et al. [2] presented numerical evidence that general sparse methods 
outperform envelope methods in both respects. However, envelope methods and related 
methods such as frontal or skyline methods continue to be the standard solution option 
in many commercial structural analysis packages. Thus, demonstrating the efficiency of 
the new spectral algorithm offers potential performance improvements in these packages 
without making substantial changes to the underlying data structures. Further, Liu [23] 
has described a generalized envelope algorithm for computing the numerical factoriza- 
tion by rows, and his results show that such a scheme can compete with general sparse 
algorithms. 

The following is an outline of the rest of this paper. In Section 2 we formulate the 
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problems associated with the minimization of envelope parameters and describe related 
problems called the 1-sum and 2-sum problems. We describe some theoretical results 
to justify the proposed new algorithm. The second Laplacian eigenvector solves a con- 
tinuous relaxation of a discrete problem related to the envelope problem, the minimum 
2-sum problem. Further, it is proved that the permutation vector computed by the 
spectral algorithm is a closest (in the 2-norm sense) permutation vector to a second 
Laplacian eigenvector. In Section 3 we discuss the spectral algorithm and its numeri- 
cal implementation. The multilevel algorithm, which uses coarsening of the underlying 
graph combined with Rayleigh Quotient iteration [3], to compute the eigenvector is de- 
scribed. Numerical results and comparisons with GPS, GK, and RCM axe presented in 
Section 4. These results indicate that the new algorithm is often considerably more effi- 
cient in reducing the storage requirements. The spectral algorithm does require greater 
execution time for computing the ordering, but the new ordering often yields greatly 
reduced factorization times for the spectrally reordered matrices. 

2. The envelope reduction problem. 

2.1. The envelope of a matrix. Let A be an n x n symmetric matrix with 
elements a tJ , whose diagonal elements are nonzero. We consider various parameters of 
the matrix A associated with its envelope. 

We denote the column indices of the nonzeros in the lower triangular part of the 
i-th row by row(i) = {j : aij ^ 0, and 1 < j < i}. For the i-th row of A we define 

(2.1) fi{A) = min{j : j (E rotn(i)}, and 

(2.2) r,(A) = i - fi(A). 

Here fi(A) is the column index of the first nonzero in the j-th row of A (by our assump- 
tion of nonzero diagonals, 1 < /, < i), and the parameter r,(A) is the row-width of the 
*-th row of A. The bandwidth of A is the maximum row-width 

bw(A ) = max{rj(A) : i = 1, . . . , n}. 

The envelope of A is the set of column indices that lie between the first nonzero 
column index and the diagonal in each row: 

Env(A) = {( i,j ) : fi(A) <j< i,and i = l,...,n}. 

We denote the size of the envelope by Esize(A) = \Env{A)\. The work in the Cholesky 
factorization of A that makes use of an envelope storage scheme can be bounded from 
above by (1/2) X^"=2 r «( r i + 3). Hence hereafter we will denote Ework(A) = X)"=i r i 35 
a measure of the work in such a factorization. We stress that this estimate is an upper 
bound on the actual work in an envelope factorization scheme. 

The values of these parameters strongly depend on the choice of an ordering of the 
rows and columns, and thus we consider how these parameters vary for a symmetrically 
permuted matrix P T AP, where P is a permutation matrix. We define Esize m i n (A), 
the minimum envelope size of A, to be the minimum size among the envelopes of all 
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permuted matrices P T AP. The quantities Ework m i„(A) and bw m i n (A) are defined in 
similar fashion. In general the minima for these three quantities will not be attained 
by the same permutation. 

The envelope parameters can also be defined with respect to the adjacency graph 
G = ( V,E ) of A. Denote nbr(v) = {v} U adj(v). In terms of the graph G and an 
ordering a of its vertices, we can define 

r(v,a ) = max{a(i;) — a(w) : w € nbr(v),at(w) < a(v)}. 


Hence we can write the envelope size and work associated with an ordering a as 

Esize(G, a) = r ( v ) = Y~1 max{a(u) — a(u>) : w 6 nbr(v), a(w) < a(v)j 
vev vzv 

Ework(G,a ) = YZ r2 ( v ) — Y"1 max{(o:(u) — a(w)) 2 : w £ nbr(v), a(w) < a(u)}. 

v£V vEV 

The goal is to choose a vertex ordering a : V {l,...,n} to minimize one of the 
parameters described above. We denote by Esize m i n (G ) (Ework mtn (G)) the minimum 
value of Esize(G,a ) (Ework min (G,a)) over all orderings a, where again (in general) 
the minima will not all be attained by the same a. We will use the definitions in terms 
of matrices throughout the rest of the paper. 

It will be helpful to consider quantities related to the envelope size and envelope 
work: the 1-sum, (?i{A), and the 2-sum, <r\{A). We write the envelope size and 1-sum, 
and the envelope work and the 2-sum in a way that shows their relationships: 


Esize(A) 

Ework(A) 

al(A) 


5Z -I 1111 * (* “ 

. =1 j€rcm/(i) 


t=l j£row(i) 


n 


E 


max (i — j) 2 , 

jerow(i) 


n 


E E.(«'-i) a - 

»= 1 j€rou/(t) 


The parameters a\ <m in{A) and cr| „(A) axe the minimum values of these parameters 
over all permuted matrices P T AP. 

It is known that minimizing the bandwidth and the 1-sum are NP-complete prob- 
lems, the former even for trees with degree bounded by three. Minimizing any of the 
other quantities considered here is likely to be intractable as well, so one has to settle 
for heuristic orderings to reduce the quantity. 

Recently it has been shown that the envelope size problem is intimately related 
to the 1-sum problem, and that the envelope work problem is related to the 2-sum 
problem [13]. Let A denote the maximum number of off-diagonal nonzeros in a row of 
A. (This is the maximum vertex degree in the adjacency graph of A.) 
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THEOREM 2.1 ([13]). Let A be a symmetric matrix . The minimum values of the 
envelope size , estimate of the envelope work in the Cholesky factorization , 1-sum, and 
2 -sum of a symmetric matrix A are related by the following inequalities : 


(2.3) 

Esizc m i n [A'j 

(2.4) 

Ework min (A) 

(2.5) 

0^2,mtn(^4) 

■ 



^ 01 ,mw(^) — A2?$tZC m j n (-A). 

— a 2,min{A) — ^Ework m i n (A). 

^ t min (^4) ^ |0’2 f m»n( J 4) • 


2.2. The Laplacian matrix and bounds on envelope parameters. The Lapla- 
cian matrix Q(G ) of an undirected graph G is the n x n matrix D — B, where D is 
the diagonal degree matrix and B is the adjacency matrix of (?. If G is the adjacency 
graph of a symmetric matrix M, then we could define the Laplacian matrix Q directly: 


~ \ 


( ~1 
0 


= i Qij 


if i ^ j and m t j 
if i ^ j and m,j = 0, 
if i = j. 


The eigenvalues of Q(G) are the Laplacian eigenvalues of G, and we list them as 
^1 < ^2 < ■ • • < ^n- An eigenvector corresponding to A* will be denoted by x k , and will 
be called a fcth eigenvector of Q. It is well-known that Q is a singular M-matrix, and 
hence its eigenvalues are nonnegative. Thus Ai = 0, and the corresponding eigenvector 
is any nonzero constant vector c. If G is connected, then Q is irreducible, and A 2 > 0. 
The smallest nonzero eigenvalues and the corresponding eigenvectors have important 
properties that make them useful in the solution of various partitioning and ordering 
problems. These properties were first investigated by Fiedler [9, 11]; more recently 
several authors have studied their application to such problems. 

Juvan and Mohar [19] have obtained bounds for bandwidth and p-sums in terms 
of Laplacian eigenvalues. They have also suggested the use of a second eigenvector to 
compute orderings to reduce bandwidth, 1-sum, and 2-sum. Helmberg, Mohar, Poljak, 
and Rendl [18] have obtained additional lower bounds on the bandwidth. The 1- and 
2-sum problems have been recently formulated as quadratic assignment problems and 
thus bounds have been obtained for the envelope size and work [13]. The following 
result describes two of the simpler bounds: 

THEOREM 2.2 ([13]). The envelope size of a symmetric matrix A can be bounded 
in terms of its second and largest Laplacian eigenvalues as 


(rc 2 - 1) < Esize min (A ) < (» 2 - !)• 

oA o 

Our estimate of the envelope work in the Cholesky factorization of A can be bounded as 


A a (A) 
12A 



- 1) < Ework min (A) < 


A n (A) 

12 


n(n 2 


- 1 ). 
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2.3. Approximate minimization of envelope work. We now offer some justi- 
fication for the spectral envelope-reduction algorithm, which computes an ordering by 
sorting the components of a second Laplacian eigenvector. The idea is to consider the 
related 2-sum problem, and then to show that a second Laplacian eigenvector x 2 solves 
a continuous relaxation of the problem. We then prove that the permutation vector 
computed by the spectral algorithm is a closest vector (in the 2-norm sense) among the 
permutation vectors to the eigenvector x 2 . 

For odd n, let V denote the set of n-vectors p whose components are permutations 
of {— (n — l)/2, ...,— 1,0, l,...,(n — l)/2). For even n, let V denote vectors that 
are permutations of {— n/2, . . . , —1, +1, . • . , n/ 2}. We denote the i-th component of a 
vector x by x t . We consider the 2-sum of a symmetric matrix A, defined with respect 
to vectors in V: 

™$£ £ (*< - = \ £ (*»■ - x i?- 

i=l j£rott;(t) 0 


A strategy to approach this hard discrete problem is to relax the condition that x must 
belong to the set of permutation vectors and instead to minimize the objective function 
over a suitable class of n-vectors. This yields an easier continuous problem; we can then 
find the permutation vector closest to the solution vector of the continuous problem, 
and consider the former as an approximate solution of the combinatorial problem. 

Note that any p£V satisfies p T u = 0, and l = p T p = (n/12)(n 2 — 1) for odd n, and 
i = (n/12)(n + l)(n + 2) for even n, where u — (1, 1, . . . 1) T . Given a vector x £ 9£ n , 
we can define a permutation vector p induced by x by the rule p, < pj if and only if 
Xi < Xj. Note that the ordering of the columns and rows is unique except when two or 
more components have the same value x,-. Hence to obtain a continuous relaxation of 
the discrete problem, we consider the set X of vectors x £ 3?" satisfying x ^ 0, xju = 0, 
and x r x = l. This is now a continuous optimization problem: 


1 . 

— min 
2 xe* 


£ (*.• - x if 

aij^O 


— min 
xex 


( > 

n 

£ dix] - 2 Y, x i x i) 

v 1=1 J<* / 

T JT / 


— min x 1 Dx — x 1 Bx = minor Qx 
xex~ ~ “ “ xex~ 


— ^ 2 ^ 2^02 


\ 2 t 


Hence a second Laplacian eigenvector x^ solves the continuous approximation to 
the 2-sum problem. Now we prove that a permutation vector p induced by x^ is a 
closest vector in V tox 2 . Earlier a similar result was obtained by Chan and Szeto [4] 
for the graph bisection problem. 

THEOREM 2.3. The vector p induced by a second Laplacian eigenvector x 2 is a 
closest (in the 2-norm) permutation vector to x 2 . In other words , 




We require the following lemma to prove the theorem. 

LEMMA 2.4. If Oi < a 2 , £>1 < b 2 are real numbers, 

r = (at - b 2 ) 2 + (a 2 - 6 X ) 2 , and s = (ai - 6 a ) 2 + (a 2 - b 2 ) 2 , 


then r > s. 

Proof. Suppose that r < s. Then 


(ai — 62 ) 2 + (02 — 61) 2 < (ai — 6j) 2 + (02 — b-f) 2 

=> a 2 (b 2 -bi) < 01(^2 — h). 

Since a\ < a 2 , it follows that b 2 < b\, which is a contradiction. ■ 

Proof of Theorem 2.3: For convenience of notation, let x = x 2 in this proof. Let 
y p m be a permutation vector such that there exists a pair of vertices u,v satisfying 
x(u) < x(v) and y(u) > y(v). Let z be the ordering such that z(u ) = y(v), z(v) = y(u ), 
and z(w) = y(w) for all other vertices. Then 

||y-x||2 2 -||z-x|| 2 2 

= (y(u) - x(u)) 2 + (y(v) - x(u)) 2 - ( y(v ) - x(u)) 2 - ( y(u ) - x(v)) 2 

> 0, 

where the last inequality follows from the previous lemma. By the swapping of com- 
ponents, we have obtained a vector z that is closer than y to the eigenvector x. By 
repeating this swapping procedure, we find that p m is a closest vector in V to the vector 
x. ■ 

Earlier Juvan and Mohar [19] had shown that p m maximized the value of the fol- 
lowing inner product over all permutation vectors p: 

|(x 2 ,p m )| > |(Sa,£)|. 

Stronger justification of the spectral algorithm for reducing the 2-sum is obtained 
in the companion paper [13] by considering a quadratic assignment formulation of the 
problem. This formulation leads to a lower bound for the 2-sum in terms of the second 
Laplacian eigenvalue, and the orthogonal matrix attaining this lower bound can be 
characterized. It can be shown that a closest permutation matrix (defined in a suitable 
sense) to this orthogonal matrix is obtained by sorting the components of a second 
Laplacian eigenvector in nondecreasing (nonincreasing) order. 

2.4. Adjacency orderings. We now consider the concept of an adjacency order- 
ing of a graph G. Let G be the adjacency graph of a matrix A, and suppose that 
the vertices of G are ordered in some ordering as {wi, . . . ,u„} (i.e., a(f j) = j), and let 
Vj = {t?i, . . . , Vj). For Y C V, define adj(Y) to be the set of vertices in V \ Y that are 
adjacent to some vertex in Y. We will say that an ordering is an adjacency ordering if 
u i+1 G adj(Vj), for j = 1, . . ., n - 1. 
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The size \adj(Vj)\ has been called the jth. frontwidth [24], and corresponds to the 
size of the j - th column of the envelope of A. Hence an alternative expression for the 
the envelope size is 


Esize(A) = \adj(Vj)\. 

3 = 1 

This expression for the envelope size shows the rationale for considering adjacency 
orderings for envelope- reduction. The idea is to locally reduce the jth frontwidth by 
choosing vj to be a vertex of low degree belonging to adj(Vj- 1 ). The Cuthill-McKee 
ordering is an adjacency ordering, but RCM is not an adjacency ordering. The GPS 
and GK algorithms attempt to number vertices in the level structures to obtain an 
adjacency ordering, as fax as is possible. 

The ordering induced by a second Laplacian eigenvector is not an adjacency order- 
ing, but comes close in the sense described below. The following theorem, proved by 
Fiedler [11], provides the necessary insight. 

THEOREM 2.5. Let G be a connected graph, and x = (xi,x 2 , . . . x n ) be a second 
Laplacian eigenvector of G. For any real p < 0, define S(p) = {v, € V : Xj > p}. Then 
the subgraph induced on S(p ) is connected. Similarly, if p > 0, then S'(p ) = {vj € V : 
Xj < p} induces a connected subgraph. 

In the notation of the theorem, let the vertices vj € V be ordered such that j < k 
if and only if xj < X*. Consider three subsets of vertices corresponding to positive, 
zero, and negative entries in the second eigenvector; i.e., define P = {vj : Xj > 0}, 
Z = {vj : Xj = 0}, and N = {vj : Xj < 0}. Let the vertices in N be numbered 
by j = 1,...,*, the vertices in Z by j = k + 1, . . ., p — 1, and the vertices in P by 
j = p, . . ., n. We have k < p. Then Theorem 2.5 implies that for j = p — 1, . . ., n, 
we have Vj + 1 G adj(Vj). A similar statement holds if we add vertices with negative 
entries in the eigenvector in decreasing order to the set PU Z. Thus the order implied 
by a second Laplacian eigenvector has the property of an adjacency ordering if vertices 
with positive components are added in increasing order to N U Z, or if vertices with 
negative components are added in decreasing order to PUZ. However, there exist simple 
examples, even trees, for which the spectral ordering is not an adjacency ordering. 

3. The Spectral algorithm for envelope reduction. Based on the theorems 
in Section 2 the following new algorithm for reducing the envelope of a sparse matrix 
can be formulated. Since the algorithm is based on properties of the spectrum of the 
Laplacian matrix L, it will be called the spectral algorithm. We assume throughout this 
section that the adjacency graph of the given matrix is connected, or that the matrix 
is irreducible. 

ALGORITHM 1. Spectral Algorithm 

1. Given the sparsity structure of a matrix M , form the Laplacian matrix L. 

2. Compute a second eigenvector x 2 of L. 

3. Sort the components of the eigenvector in nondecreasing order, and reorder the 
matrix M using the corresponding permutation vector. Also sort the compo- 
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nents in nonincreasing order, and compute the corresponding reordering of the 
matrix M. Choose the permutation that leads to the smaller envelope size. 

The implementation of steps 1 and 3 are relatively straightforward. The formation 
of the Laplacian matrix requires the computation of the degree of the nodes x,. Step 
3 is a simple sort of the entries of x 2 , and recording the resulting permutation of 
indices. This can be done quickly by any efficient sorting algorithm such as quicksort. 
Computationally the difficult part is step 2. 

The standard algorithm for computing a few eigenvalues and eigenvectors of large 
sparse symmetric matrices is the Lanczos algorithm. Since the Lanczos algorithm is 
discussed extensively in the textbook literature [16, 26], we do not include a detailed 
description of the standard algorithm here. Recently, we have developed a much more 
efficient multilevel method for finding a second eigenvector [3]. The multilevel method 
requires three elements in addition to the Lanczos algorithm: 

• Contraction: Construct a series of smaller graphs that in some sense retain 
the global structure of the original large graph. 

• Interpolation: Given a second eigenvector of a contracted graph, interpolate 
this vector to the next larger graph in a way that provides a good approximation 
to an eigenvector of the larger graph. 

• Refinement: Given an approximate eigenvector for a graph, compute a more 
accurate vector efficiently. 

Graph contraction is accomplished by first finding a maximal independent set of ver- 
tices, which are to be the vertices of the contracted graph. The edges of the contracted 
graph are determined by growing domains from the selected vertices in a breadth-first 
manner, adding an edge to the contracted graph when two domains intersect. A series 
of smaller contracted graphs is constructed until the size of the vertex set is less than 
some number (typically 100). The Lanczos algorithm can then be used to find the 
eigenvector of the smallest graph very quickly. This eigenvector is then interpolated 
to a vector corresponding to the next larger graph. This interpolated vector yields a 
very good approximation to the eigenvector of the larger graph. The approximation 
is then refined using the Rayleigh Quotient Iteration algorithm, which, because of its 
cubic convergence, usually requires only one or perhaps two iterations to obtain an 
acceptable result. This process of interpolation and refinement is continued until the 
eigenvector of the original graph is determined. 

4. Numerical results. This section shows numerical results for the envelope sizes 
and bandwidths obtained from the spectral, RCM, GPS, and GK algorithms for three 
sets of matrices. The first set, shown in Table 4.1, includes matrices for structural 
analysis applications from the Boeing-Harwell data set. The next set, shown in Ta- 
ble 4.2, consists of miscellaneous matrices from the Boeing-Harwell collection. Finally, 
the third set, shown in Table 4.3, is a selection of matrices from structural analysis used 
at NASA. The computations were performed on a Silicon Graphics workstation with a 
33 MHZ IP7 processor. 

The spectral algorithm finds the reordering with the smallest envelope in 14 out of 
18 cases (as shown in the “Rank” column of the tables). In those cases in which the re- 

9 




Table 4.1 

Results (Boeing- Harwell — Structural Analysis) 


Title 

(equations) 

(nonzeros) 

Envelope 

Bandwidth 

Run time 
(sec.) 

Algorithm 

Rank 

BCSSTK13 

64,486 

455 

3.92 

SPECTRAL 

4 

(2,003) 

58,542 

223 

.64 

GK 

3 

(11,973) 

57,501 

145 

.57 

GPS 

2 


56,299 

198 

.08 

RCM 

1 

BCSSTK29 

3,067,004 

882 

31.95 

SPECTRAL 

1 

(13,992) 

6,948,091 

1,505 

9.53 

GK 

2 

(316,740) 

7,040,998 

869 

5.29 

GPS 

3 


7,374,140 

914 

2.37 

RCM 

4 

BCSSTK30 

9,135,742 

4,769 

78.18 

SPECTRAL 

1 

(28,924) 

15,686,968 

16,947 

78.10 

GK 

2 

(1,036,208) 

23,242,990 

2,515 

61.65 

GPS 

3 


23,242,990 

2,512 

6.32 

RCM 

4 

BCSSTK31 

19,574,992 

4,763 

55.06 

SPECTRAL 

1 

(35,588) 

22,330,987 

1,880 

22.05 

GK 

2 

(608,502) 

23,416,579 

1,104 

9.12 

GPS 

3 


23,641,124 

1,176 

4.69 

RCM 

4 


27,614,531 

13,792 

92.09 1 

SPECTRAL 

nsy 


49,457,764 

3,761 

102.44 

GK 

11 


50,067,390 

2,339 

79.48 

GPS 



52,170,122 

2,390 

7.83 

RCM 

■1 

MswaaiMreTil 

3,788,702 

1,199 

31.01 

SPECTRAL 

3 

(8,738) 

3,571,395 

932 

5.20 

GK 

1 


3,717,032 

519 

3.22 

GPS 

2 


3,799,285 

749 

1.82 

RCM 

4 
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Table 4.2 

Results (Boeing-Harwell — Miscellaneous) 


Title 

(equations) 

(nonzeros) 

Envelope 

Bandwidth 

Run time 
(sec.) 

Algorithm 

Rank 

CAN1072 

55,228 

301 

.51 

SPECTRAL 

2 

(1,072) 

48,538 

234 

.20 

GK 

1 


74,067 

159 

.13 

GPS 

4 


56,361 

175 

.05 

RCM 

3 

POW9 

29,149 

264 

.45 

SPECTRAL 

1 

1 . 

64,788 

201 

.14 

GK 

2 


69,446 

116 

.10 

GPS 

3 

■Him 

79,260 

133 

.05 

RCM 

4 

BLKHOLE 

120,767 

426 

.56 

SPECTRAL 

1 

(2,132) 

169,219 

134 

.17 

GK 

2 

(8,502) 

173,243 

■ 

.12 

GPS 

4 


171,437 



RCM 

3 


93,907 

142 

.78 

SPECTRAL 

1 


96,591 

92 

.28 

GK 

2 

(13,853) 

101,769 

65 

.19 

GPS 

3 


102,983 

69 

.11 

RCM 

4 

SSTMODEL 

86,635 


2.21 

SPECTRAL 

1 

(3,345) 



.28 

GK 

2 

(13,047) 

110,936 

83 

.17 

GPS 

4 


105,421 

88 

.10 

RCM 

3 
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Table 4.3 
Results (NASA) 


Title 

(equations) 

(nonzeros) 

Envelope 

Bandwidth 

Run time 
(sec) 

Algorithm 

Rank 

BARTH4 

345,623 

593 

1.60 

SPECTRAL 

1 

(6,019) 

658,181 

280 

.54 

GK 

2 

(23,492) 

669,239 

213 

.33 

GPS 

3 


725,950 

215 

.21 

RCM 

4 

SHUTTLE 

566,496 

631 

2.59 

SPECTRAL 

3 

(9,205) 

531,420 

92 

1.12 

GK 

1 

(45,966) 

531,422 

92 

.93 

GPS 

2 


567,887 

150 

.32 

RCM 

4 

SKIRT 

688,924 

1,021 

5.14 

SPECTRAL 

1 

(12,598) 

1,013,423 

425 

3.20 

GK 

2 

(104,559) 

1,039,544 

309 

2.46 

GPS 

3 


1,068,993 

314 

.82 

RCM 

4 

PWT 

5,101,527 

1,627 

13.62 

SPECTRAL 

1 


5,520,603 

450 

29.65 

GK 

2 

■IE* ft' 

5,638,855 

340 

28.27 

GPS 

4 

mmm 

5,652,184 

340 

1.67 

RCM 

3 

BODY 

6,706,747 

2,496 

26.60 

SPECTRAL 

1 

(45,087) 

10,526,446 

1,081 

13.60 

GK 

2 

(208,821) 

10,658,164 

667 

8.42 

GPS 

3 


11,470,411 

756 

2.23 

RCM 

4 

FLAP 

10,471,456 

1,784 

45.90 

SPECTRAL 

1 

(51,537) 

12,367,171 

1,019 

24.96 

GK 

3 

(531,157) 

12,339,642 

743 

19.08 

GPS 

2 


12,598,705 

874 

4.19 

RCM 

4 

IN3C 

425,232,466 

9,504 

117.83 

SPECTRAL 

1 

(262,620) 

519,316,395 

3,780 

56.97 

GK 

2 

(1,026,888) 

526,302,263 

2,473 

26.28 

GPS 

3 


581,700,745 

2,746 

12.88 

RCM 

4 
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suit of the spectral algorithm is not the best (i.e., BCSSTK13, BKSSTK33, SHUTTLE, 
and CAN 1072), it is still fairly close to the best result. In several cases, however, the 
spectral algorithm finds a reordering with an envelope substantially smaller than any 
of the other algorithms, sometimes by a factor of more than two. Note also that the 
spectral algorithm clearly outperforms the others on the larger problems in the Tables. 
The run time of the spectral algorithm is usually, but not always, greater than that 
of the other algorithms. We expect the differences in runtimes between the ordering 
algorithms to be smaller on computers with vector-processing capabilities, such as the 
Crays. 

The GPS, GK, and RCM algorithms, which are all closely related, use local search 
(breadth-first search) from a pseudo-peripheral vertex to generate a long rooted level 
structure. The RCM algorithm then numbers the vertices by increasing level values, 
where the vertices in each level axe numbered in nondecreasing order of their degrees. 
The final RCM ordering is obtained by reversing the ordering thus obtained. The GPS 
and GK algorithms use more sophisticated techniques to create a more general level 
structure by combining the information from two rooted level structures obtained from 
the endpoints of a pseudo-diameter in the RCM algorithm. They also use more refined 
numbering techniques to reduce the size of the envelope and the bandwidth. This is 
the reason why the latter two algorithms require more time than the RCM algorithm. 

Generally the GPS algorithm yields a lower bandwidth while the GK algorithm 
yields a lower envelope size [14, 21]. Our results are in agreement with this conclusion. 
It should be pointed out that n = 2680 was the largest order of the problems considered 
in earlier work, and that the results reported here are for much larger problems. 

In contrast to the above algorithms, the spectral algorithm relies on the global 
information in the components of a second Laplacian eigenvector. The results show 
that the bandwidths of the spectral reorderings are often much greater than those of 
the other reorderings, even when the spectral envelopes are much smaller. This can be 
seen in Figures 4.1 through 4.5, which show the sparse matrix structure of the original 
BARTH4 matrix and of the four reorderings considered here. A black dot indicates a 
nonzero element. The GK, GPS, find RCM reorderings all look very similar, whereas 
the SPECTRAL reordering has a quite different appearance. 

Table 4.4 

Factorization times 


Title 

Envelope 

Factor time 
(sec) 

Algorithm 

BCSSTK29 

3,067,004 


SPECTRAL 


7,374,140 


RCM 

BCSSTK33 

3,788,702 

670 



3,799,285 

685 


BARTH4 

345,623 




725,950 




Juvan and Mohar [19] had suggested the use of the spectral ordering for reducing 
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the bandwidth (and p-sums), but our results show that the GPS algorithm is much 
more effective than the spectral algorithm in reducing the bandwidth. A possibility is 
to make limited use of a local reordering strategy based on the adjacency structure to 
improve the envelope parameters obtained from the spectral method. Such reordering 
strategies will be considered elsewhere since the evaluation of the various possibilities 
will require much effort. 

Finally we list in Table 4.4 the factorization times for a few matrices, reordered 
with both the spectral algorithm and with RCM. These times are for the envelope 
factorization routine from SPARSPAK, and are measured again on a SGI worksta- 
tion. We selected one example where the spectral algorithm is comparable in storage 
requirements to RCM (BCSSTK33), and two examples where the spectral algorithm 
yields considerably lower storage memory requirements. The results demonstrate the 
quadratic behavior of the factorization time as a function of the envelope size. There- 
fore we conclude that spectral reordering not only reduces the memory requirements, 
but also improves execution times. 
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